Notebook source code: notebooks/03_methods_estimate_manifold_topology.ipynb

Estimate Neural Topology

Set Up + Imports

import setup

setup.main()
%load_ext autoreload
%autoreload 2
%load_ext jupyter_black

import matplotlib.pyplot as plt
import numpy as np
import plotly.graph_objects as go
from gtda.diagrams import PairwiseDistance
from gtda.homology import WeakAlphaPersistence
from gtda.plotting import plot_diagram, plot_heatmap
from plotly.subplots import make_subplots

# import neurometry.curvature.datasets.experimental as experimental
import neurometry.curvature.datasets.gridcells as gridcells
import neurometry.topology.persistent_homology as persistent_homology
import neurometry.curvature.viz as viz
from neurometry.curvature.datasets.synthetic import (
    load_s2_synthetic,
    load_t2_synthetic,
)

Working directory:  /home/facosta/neurometry/neurometry
Directory added to path:  /home/facosta/neurometry
Directory added to path:  /home/facosta/neurometry/neurometry
The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
The jupyter_black extension is already loaded. To reload it, use:
  %reload_ext jupyter_black

Persistence homology for synthetic sphere, torus point clouds

Generate point clouds

dim = 3

sphere_point_cloud, sphere_labels = load_s2_synthetic(
    "random",
    n_times=400,
    radius=1,
    geodesic_distortion_amp=0,
    embedding_dim=dim,
    noise_var=0.001,
)
sphere_point_cloud = np.array(sphere_point_cloud)

torus_point_cloud, torus_labels = load_t2_synthetic(
    "random",
    n_times=400,
    major_radius=2,
    minor_radius=1,
    geodesic_distortion_amp=0,
    embedding_dim=dim,
    noise_var=0.0001,
)

torus_point_cloud = np.array(torus_point_cloud)

Plot point clouds

s2_x, s2_y, s2_z = sphere_point_cloud.T

t2_x, t2_y, t2_z = torus_point_cloud.T

fig = make_subplots(
    rows=1, cols=2, specs=[[{"type": "scatter3d"}, {"type": "scatter3d"}]]
)


fig.add_trace(
    go.Scatter3d(
        x=s2_x, y=s2_y, z=s2_z, mode="markers", marker=dict(color="red", size=3)
    ),
    row=1,
    col=1,
)

fig.add_trace(
    go.Scatter3d(
        x=t2_x, y=t2_y, z=t2_z, mode="markers", marker=dict(color="orange", size=3)
    ),
    row=1,
    col=2,
)

fig.update_layout(
    width=1000,
    height=500,
    title_text="3D Scatter Plots of Sphere and Torus Point Clouds with Smaller Points",
)

fig.show()

sphere_point_cloud.shape

(400, 3)

# Compute the 0, 1 and 2-dimensional homology of the point clouds
homology_dimensions = (0, 1, 2)
WA = WeakAlphaPersistence(homology_dimensions=homology_dimensions)

# Compute the persistence diagrams
diagrams = WA.fit_transform([sphere_point_cloud, torus_point_cloud])
print(diagrams.shape)

(2, 726, 3)

fig_sphere = plot_diagram(
    diagrams[0],
    homology_dimensions=(0, 1, 2),
    plotly_params={"title": "Sphere"},
)
fig_sphere.update_layout(title="Sphere")

Data type cannot be displayed: application/vnd.plotly.v1+json

fig_torus = plot_diagram(
    diagrams[1],
    homology_dimensions=(0, 1, 2),
    plotly_params={"title": "Torus"},
)

fig_torus.update_layout(title="Torus")

Data type cannot be displayed: application/vnd.plotly.v1+json

Persistence diagrams for synthetic 2-torus

diagrams = persistent_homology.compute_persistence_diagrams(
    torus_point_cloud, maxdim=2, n_threads=-1
)

'compute_persistence_diagrams' executed in 41.9059s

viz.plot_persistence_diagrams(diagrams, density=True)

Compute Bottleneck Distance between Persistence Diagrams

Create torus point clouds with varying levels of noise

noise_var = np.logspace(-5, 0, 15)

tori = []

for noise in noise_var:
    torus, _ = load_t2_synthetic(
        None,
        n_times=1800,
        major_radius=2,
        minor_radius=1,
        geodesic_distortion_amp=0,
        embedding_dim=dim,
        noise_var=noise,
    )
    tori.append(torus)

Compute persistence diagrams for tori and compute pairwise bottleneck distance

# Compute the 0, 1 and 2-dimensional homology of the torus
homology_dimensions = (0, 1, 2)
WA = WeakAlphaPersistence(homology_dimensions=homology_dimensions)

# Compute the diagrams for the tori
diagrams = WA.fit_transform(tori)
print(diagrams.shape)

# Compute the bottleneck distance between the diagrams
PD = PairwiseDistance(metric="bottleneck", n_jobs=-1)

distance = PD.fit_transform(diagrams)
print(distance.shape)

(15, 3747, 3)
(15, 15)

plt.scatter(noise_var, distance[0, :], label="0D")
# make the plot look nice
plt.xscale("log")
plt.yscale("log")
plt.xlabel("Noise Variance")
plt.ylabel("Bottleneck Distance to Original Torus")

plt.grid()

plot_heatmap(distance, colorscale="blues")

Data type cannot be displayed: application/vnd.plotly.v1+json

# diagram_0 = persistent_homology.compute_persistence_diagrams(
#     torus_0, maxdim=2, n_threads=-1
# )

# diagram_1 = persistent_homology.compute_persistence_diagrams(
#     torus_1, maxdim=2, n_threads=-1
# )

'compute_persistence_diagrams' executed in 46.4076s
'compute_persistence_diagrams' executed in 31.8140s

Persistence homology for place cell data

Load place cell data

From Ravikrishnan P Jayakumar, Manu S Madhav, Francesco Savelli, Hugh T Blair, Noah J Cowan, and James J Knierim. Recalibration of path integration in hippocampal place cells. Nature, 566(7745):533–537, 2019.

expt_id = 34
timestep = int(1e6)

dataset, labels = experimental.load_neural_activity(
    expt_id=expt_id, timestep_microsec=timestep
)
dataset = dataset[labels["velocities"] > 5]
labels = labels[labels["velocities"] > 5]
dataset = np.log(dataset.astype(np.float32) + 1)
dataset.shape

INFO: # - Found file at /home/facosta/neurometry/neurometry/data/binned/expt34_times_timestep1000000.txt! Loading...
INFO: # - Found file at /home/facosta/neurometry/neurometry/data/binned/expt34_place_cells_timestep1000000.npy! Loading...
INFO: # - Found file at /home/facosta/neurometry/neurometry/data/binned/expt34_labels_timestep1000000.txt! Loading...

(934, 40)

Persistence diagrams for place cell data

place_cell_diagrams = compute_persistence_diagrams(dataset, maxdim=2, n_threads=-1)
plot_persistence_diagrams(place_cell_diagrams)

Function 'compute_persistence_diagrams' executed in 1.6295s

Synthetic Grid cell data

Generate synthetic grid cell data + compute persistence diagrams

Orientation variability = 0

grid_scale = 1
arena_dims = np.array([4, 4])
n_cells = 256
grid_orientation_mean = 0
grid_orientation_std = 0

field_width = 0.05
resolution = 50

neural_activity, _ = gridcells.load_grid_cells_synthetic(
    grid_scale,
    arena_dims,
    n_cells,
    grid_orientation_mean,
    grid_orientation_std,
    field_width,
    resolution,
)
print("shape of neural activity matrix: " + str(neural_activity.shape))

shape of neural activity matrix: (2500, 256)

diagrams = compute_persistence_diagrams(neural_activity, maxdim=2, n_threads=-1)
plot_persistence_diagrams(diagrams)

Function 'compute_persistence_diagrams' executed in 4809.5422s

grid_scale = 1
arena_dims = np.array([4, 4])
n_cells = 256
grid_orientation_mean = 0
grid_orientation_std = 3

field_width = 0.05
resolution = 50

neural_activity, _ = gridcells.load_grid_cells_synthetic(
    grid_scale,
    arena_dims,
    n_cells,
    grid_orientation_mean,
    grid_orientation_std,
    field_width,
    resolution,
)
print("shape of neural activity matrix: " + str(neural_activity.shape))

shape of neural activity matrix: (2500, 256)

Shuffle Data and plot diagrams

import os

os.environ["GEOMSTATS_BACKEND"] = "pytorch"
import geomstats.backend as gs

import neurometry.datasets.synthetic as synthetic

task_points = synthetic.hypersphere(1, 1000)
noisy_points, manifold_points = synthetic.synthetic_neural_manifold(
    points=task_points,
    encoding_dim=3,
    nonlinearity="sigmoid",
    scales=gs.array([1, 1, 1]),
    poisson_multiplier=100,
)

noise level: 0.71%

!pip install dreimac

Collecting dreimac
  Downloading dreimac-0.3.1-py3-none-any.whl.metadata (6.2 kB)
Requirement already satisfied: matplotlib>=3.6 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from dreimac) (3.8.4)
Requirement already satisfied: numba>=0.56 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from dreimac) (0.59.1)
Requirement already satisfied: numpy>=1.23 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from dreimac) (1.26.4)
Requirement already satisfied: persim>=0.3 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from dreimac) (0.3.5)
Requirement already satisfied: ripser>=0.6 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from dreimac) (0.6.8)
Requirement already satisfied: scipy>=1.10 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from dreimac) (1.13.0)
Requirement already satisfied: contourpy>=1.0.1 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from matplotlib>=3.6->dreimac) (1.2.1)
Requirement already satisfied: cycler>=0.10 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from matplotlib>=3.6->dreimac) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from matplotlib>=3.6->dreimac) (4.51.0)
Requirement already satisfied: kiwisolver>=1.3.1 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from matplotlib>=3.6->dreimac) (1.4.5)
Requirement already satisfied: packaging>=20.0 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from matplotlib>=3.6->dreimac) (24.0)
Requirement already satisfied: pillow>=8 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from matplotlib>=3.6->dreimac) (10.3.0)
Requirement already satisfied: pyparsing>=2.3.1 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from matplotlib>=3.6->dreimac) (3.1.2)
Requirement already satisfied: python-dateutil>=2.7 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from matplotlib>=3.6->dreimac) (2.9.0.post0)
Requirement already satisfied: llvmlite<0.43,>=0.42.0dev0 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from numba>=0.56->dreimac) (0.42.0)
Requirement already satisfied: scikit-learn in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from persim>=0.3->dreimac) (1.4.2)
Requirement already satisfied: hopcroftkarp in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from persim>=0.3->dreimac) (1.2.5)
Requirement already satisfied: deprecated in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from persim>=0.3->dreimac) (1.2.14)
Requirement already satisfied: joblib in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from persim>=0.3->dreimac) (1.4.0)
Requirement already satisfied: Cython in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from ripser>=0.6->dreimac) (0.29.37)
Requirement already satisfied: six>=1.5 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from python-dateutil>=2.7->matplotlib>=3.6->dreimac) (1.16.0)
Requirement already satisfied: wrapt<2,>=1.10 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from deprecated->persim>=0.3->dreimac) (1.16.0)
Requirement already satisfied: threadpoolctl>=2.0.0 in /home/facosta/miniconda3/envs/neurometry/lib/python3.11/site-packages (from scikit-learn->persim>=0.3->dreimac) (3.4.0)
Downloading dreimac-0.3.1-py3-none-any.whl (36 kB)
Installing collected packages: dreimac
Successfully installed dreimac-0.3.1

from dreimac import CircularCoords
from persim import plot_diagrams

# prepare plot with 4 subplots
f, (a0, a1, a2, a3) = plt.subplots(1, 4, width_ratios=[1, 1, 1, 0.2], figsize=(14, 3))


# plot the persistence diagram, showing a single prominent class
cc = CircularCoords(X, n_landmarks=200)
plot_diagrams(cc._dgms, title="Persistence diagram", ax=a1)

# plot the data colored by the circle-valued map constructed by DREiMac
circular_coordinates = cc.get_coordinates()
a2.scatter(X[:, 0], X[:, 1], c=circular_coordinates, s=10, cmap="viridis")
a2.set_title("Input colored by circular coordinate")
a2.axis("off")
a2.set_aspect("equal")

# plot colorbar
img = a3.imshow([[0, 1]], cmap="viridis")
a3.set_visible(False)
cb = plt.colorbar(mappable=img, ticks=[0, 0.5, 1])
_ = cb.ax.set_yticklabels(["0", "180", "360"])

# plot in 3d
fig = go.Figure()
fig.add_trace(
    go.Scatter3d(
        x=task_points[:, 0],
        y=task_points[:, 1],
        z=noisy_points[:, 2],
        mode="markers",
        marker=dict(size=3),
    )
)
fig.show()

X = noisy_points


def shuffle_entries(data):
    # Shuffle each row's entries independently
    shuffled_data = np.apply_along_axis(np.random.permutation, 1, data)
    return shuffled_data


X_shuffled_1 = shuffle_entries(X)
X_shuffled_2 = shuffle_entries(X)

fig = go.Figure()
fig.add_trace(
    go.Scatter3d(
        x=X_shuffled_1[:, 0],
        y=X_shuffled_1[:, 1],
        z=X_shuffled_1[:, 2],
        mode="markers",
        marker=dict(size=3),
    )
)

fig.add_trace(
    go.Scatter3d(
        x=X_shuffled_2[:, 0],
        y=X_shuffled_2[:, 1],
        z=X_shuffled_2[:, 2],
        mode="markers",
        marker=dict(size=3),
    )
)

rng = np.random.default_rng()

# n_permutations = 10

X_shuff_1 = rng.permutation(X, axis=0)
X_shuff_2 = rng.permutation(X, axis=0)

fig = go.Figure()
fig.add_trace(
    go.Scatter3d(
        x=X_shuff_1[:, 0],
        y=X_shuff_1[:, 1],
        z=X_shuff_1[:, 2],
        mode="markers",
        marker=dict(size=3),
    )
)

fig.add_trace(
    go.Scatter3d(
        x=X_shuff_2[:, 0],
        y=X_shuff_2[:, 1],
        z=X_shuff_2[:, 2],
        mode="markers",
        marker=dict(size=3),
    )
)

from neurometry.datasets.synthetic import hypertorus, synthetic_neural_manifold

num_points = 1000
intrinsic_dim = 2
encoding_dim = 1000

torus_points = hypertorus(intrinsic_dim, num_points)

noisy, manifold = synthetic_neural_manifold(torus_points, encoding_dim, "linear")

WARNING! Poisson spikes not generated: mean must be non-negative
noise level: 7.07%

print(manifold.shape)

from neurometry.topology.persistent_homology import compute_persistence_diagrams

diagrams = compute_persistence_diagrams([manifold])

plot_diagram(diagrams[0], homology_dimensions=(0, 1, 2))

torch.Size([1000, 1000])

transposed_manifold = manifold.T
print(transposed_manifold.shape)
transposed_diagrams = compute_persistence_diagrams([transposed_manifold])

plot_diagram(transposed_diagrams[0], homology_dimensions=(0, 1, 2))

torch.Size([1000, 1000])

from sklearn.decomposition import PCA

pca = PCA(n_components=10)
pca_manifold = pca.fit_transform(manifold)

pca_diagrams = compute_persistence_diagrams([pca_manifold])

plot_diagram(pca_diagrams[0], homology_dimensions=(0, 1, 2))

pca_transposed_manifold = pca.fit_transform(transposed_manifold)

pca_transposed_diagrams = compute_persistence_diagrams([pca_transposed_manifold])

plot_diagram(pca_transposed_diagrams[0], homology_dimensions=(0, 1, 2))